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Abstract 

The switching probability of a single-domain ferromagnet under spin-current excitation is evalu- 
ated using the Fokker-Planck equation(FPE). In the case of uniaxial anisotropy, the FPE reduces to 
an ordinary differential equation in which the lowest eigenvalue Ai determines the slowest switching 
events. We have calculated Ai by using both analytical and numerical methods. It is found that the 
previous model based on thermally distributed initial magnetization states [lj] can be accurately 
justified in some useful limiting conditions. 
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I. INTRODUCTION 



Fast and reliable nanosecond level writing is essential for spin-torque-induced switching 
in memory and recording technologies. One intrinsic source for write threshold distribution 
is thermal fluctuation, two aspects of which could affect reversal time r s . First, the initial 
position of the magnetic moment is thermally distributed at the time the reversal field or 
current is applied, causing a variation in switching time. Secondly, during reversal, thermal 
fluctuation would modify the orbit, causing additional fluctuation for t s even for identical 
initial conditions. A reversal time distribution due to a thermalized initial condition was 
recently estimated [1|]. Here we present a Fokker-Planck formulation for the full problem 
including both thermalized initial condition and reversal orbit with estimates for the reversal 
time and its distribution. 

Fokker-Planck equations have been successfully applied to macro-spin systems with a 
static thermal distribution for the extraction of a small thermal escape probability. For ex- 
ample, the distribution for a thermally activated reversal has been analytically and numeri- 
cally obtained when the applied magnetic field is smaller than the coercive field 
There the switching is through thermally assisted transition from one metastable energy 
minimum to equilibrium. The ensemble-averaged t s depends exponentially on the ratio 
between the reversal energy barrier and thermal energy k^T . These subthreshold results 
have recently been extended to include the effect of spin-torque, which modifies the effective 
temperature [7, 8, 9]. 

For fast, dynamically driven switching, the applied spin-torque is above the zero- 
temperature reversal threshold. For finite temperature, this leads to a time-dependent 
Fokker-Planck equation in which the time-evolution of the probability represents the ther- 
malized reversal of the macro-spin. Here we solve this problem in the colinear geometry 
with a uniaxial anisotropy energy landscape and the presence of a spin torque. 



II. FOKKER-PLANCK EQUATION 

Consider an ensemble time-dependent magnetization probability density P(n m ,t), where 
n m is the unit vector describing the direction of the macro-spin moment m; in spherical 
coordinates n m is characterized by (9,4>). Before one turns on the magnetic field and the 
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current, the probability density takes its equilibrium value. For a uniaxial anisotropy-only 
situation, 

P(n m ,0) = P exp(-£sin 2 #) (1) 

for < 9 < 7r/2 and P(n m ,0) = for tt/2 < 6 < tt, where Po is the normalization factor 
(J P sin 6d0 = 1), Po ~ 2£, if we consider £ = KV/ksT ^> 1 where if is the anisotropy 
constant and V is the volume of nanomagnet. To determine P(n m , t) after one turns on the 
field and the current at t = 0, we need to solve the Fokker-Planck equation given below 



aT + v ' J 



DV 2 P = 



(2) 



where J is the probability current and D = ayksT/m is the diffusion constant, a is the 
damping parameter, 7 is the gyromagnetic coefficient, and m = |m|. The probability current 
is 



P- 



dt 



-7Pn m x [H e + an m x (H e - H s )] 



(3) 



where we have used the Landau-Lifshitz-Gilbert equation including the spin torque term 
H s = (Ip)(h/2q)(l/am)n s (I is current density, p is the spin polarization coefficient, q is 
the electron charge, and n s is a unit vector in the direction of the magnetization of the 
pinned layer). To make the solution mathematically tractable, we consider the case where 
the magnetic field H and n s are parallel to the anisotropy axis, i.e., H e = (H + HkCos9)e z 
and n s = e z where = 2K/M S . With these simplifications, Eq. (3) reduces to 
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where x = cos 8, and we have defined the effective potential U e = (H s — H)mx — (Hk/2)mx 2 . 
Equation (4) can be solved via the method of separation of variables. Let's assume P(x, t) = 
f(t)u(x) and it is easy to see that f(t) = exp(— Xt) and u(x) satisfies 
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The above equation can be converted into a Sturm-Liouville equation, 
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where F(x) = e /3Ue ^u(x) and c = \m/ a^ksT . Equation (6) is the same equation introduced 
by Brown except that U e includes the spin-current term H s . The original FPE (4) is 



now reduced to the standard eigenvalve problem. Namely, we determine eigen-function 
F(x) = F n (x) and eigenvalue A = X n from Eq. (6) for n = 0, 1, 2, and the general solution 
of Eq. (4) is 

P( x , t)=J2 A n e^ u ^F n (x)e~ Xnt (7) 

n=0 

where the coefficients in Eq. (7) are determined by the initial condition Eq. (1): A n = 
f dxP(x, 0)F n (x), where we have chosen the eigenfunction F n {x) to be weighted normal- 
ized: $l x dxe- pu *(*)Fl(x) = \. 

The remaining task is to solve eigenvalue A n (or c n ) and eigenfunction F n . Obviously, 
Ao = is the lowest eigenvalue with the corresponding eigenfunction Fq being a constant. 
Thus, u (x) oc exp(— (3U e ) describes the thermal equilibrium state. The smallest nonzero 
eigenvalue Ai then determines the slowest decaying or switching speed. 

The early theoretical attempts toward solving Eq. (6) 0, [3, y], 5, (jjliad been focused on 
the cases where there are two potential wells separated by an energy barrier and H < H^. 
These solutions do not apply to the fast switching case where H > H^- In this case, there 
is no energy barrier and the only stable solution is at x = — 1 or 9 = n. In the following, 
we develop two methods to solve the Eq. (6): one is the approximate minimization which 
provides an upper bound for Ai and the other is to numerically determine the coefficients of 
polynomial series when Eq. (6) is expanded into polynomials. 

The variational method can be obtained by multiplying F(x) to Eq.(6) and then by 
integrating from -1 to 1, 



1 , r,d 

dxF 

-l 



dx 



'l-x 2 )e 
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dF 
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l 

< / e~ l3Ue F 2 dx = 0. 



Performing integration by parts of the first term, we find 



c = V/H 



(8) 



where 



V = f\dxe-^(l-x 2 ) 



(9) 



H = jt 1 dxQ-^F 2 {x) 

Since both T> and 7i are definitely positive, the variational principle is to choose a trial 
function F(x) such that c is minimized, and the non-zero minimum c represents the upper 
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bound of the lowest eigenvalue c\. Additional constriction of F(x) is its orthogonality with 
F which is a constant, i.e., 

J e-< 3Ue F(x)dx = Q. (10) 

Since exp(— f3U e ) is sharply peaked at x — — 1 in our case, we change the variable as x = y — 1 
and thus f3U e = rjy — £y 2 + const, where we define rj = -^(H s — H + H k ) and £ = -j^f\H k . 
We consider a simple trial function F(y) = 1 + b±y + b 2 y 2 where b\ and b 2 are variational 
parameters. By placing F(y) into Eq. (9) and (10), we have 

V = f*y(2 - y)e^y\b 1 + 2b 2 yfdy 

n = J 2 e-w+^l + b lV + b 2 y 2 ) 2 dy (11) 

k ^e- r 'y^y 2 (l + b iy + b 2 y 2 )dy = 

After the completion of these three integrations, the eigenvalue c\ is obtained by minimizing 
V/Ti with respect to b 2 or b\. The above integrations can be analytically obtained if we we 
notice that y is limited to 1/n due to exponential e~ vy and thus the term e^ y2 is small and 
can be expanded as: e^ 2 = 1 + £y 2 + \i 2 y i - By using J Q 2 y n e~ r,y dy ~ n\/n n+1 , for example, 
the third equation in Eq. (11) is 

/ 2 24£ 360£ 2 \ , 

Note that the leading orders of b\ and b 2 are 77 and rf respectively, and the above equation 
is valid up to i]~ 3 . If we consider the fast switching case, i.e. n ^> £, or H s — H ^> the 
equation (12) can be further approximated as 

61 = -2&2??- 1 - 77 (13) 

We can similarly calculate 7i and P. Substituting b\ by Eq. (13), the expression of V/H 
contains only b 2 . Then, we find b 2 by minimizing T>/7i 

d (V 
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The results are: 

h„ = 

2 



b 2 = bi = -2 V 



n = \, V = 4, Cl = % = A V 
Rewriting in terms of original parameters, we find the slowest decaying rate is 

Ai = 4« 7 (H s -H) = 4(/// c - IK" 1 (15) 

where To" 1 = ^f/ c = ai {H + H k ). 

Although the variational method is able to analytically obtain the upper bound for the 
rate of the relaxation, its accuracy is unknown. A more rigorous method is to expand F(x) 
via Taylor series, i.e., F(x) = J2a n x n (n = 0, 1,2...). By substituting it into Eq. (6), we find 
the recurrence formula for coefficient a n 

(n + 2)(n + l)a n+2 - CO + l)«Wi 
+ [c n - n(n - 1) + n(2£ - 2)}a n 
+((n - l)a n _i - 2^{n - 2)a„_ 2 = (16) 

where ( = r\ — 2£ = ^f{H s — H). In the case of ( = £ = 0, the above equation reduces 
to a two-term recursion formula, and the solutions are the Legendre polynomials with c n = 
n(n + 1). This is the well-known solution for the thermal decay of a free particle 2|. In the 
case of ( = 0, Eq. (16) reduces to a three-term recursion formula, and the continued fraction 
method can be used to obtain the smallest nonzero eigenvalue c\ = 2 — |£ + 0(^ 2 )[3|. 

In our case where both £ and £ are non-zero, above methods can not be used directly. 
One way to solve the problem is to numerically calculate Eq. (16) by simply keeping enough 
terms (i.e. up to a large number iV = 300) to ensure the convergence of c\. In Fig. l,we show 
that the relaxation rates as a function of temperature (Hk is fixed) for different J// c (> 1). 
In Fig. 2, we show that the relaxation rate is linearly dependent of J/J c — 1 when I/I c > 1.5. 
By comparing with Eq. (15) derived from the variational approach, we find that to a good 
approximation the relaxation rate is a factor of two smaller, i.e., 

X 1 = 2(I/I c -l)r 1 . (17) 
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III. COMPARISON WITH SUN'S MODEL 

Sun's model [lj assumes that thermally distributed initial magnetization states determine 
the distribution of switching time in case of / > I c . In this model, the switching time at 
zero temperature is estimated as r = r (J// c — l) -1 ln(7r/2#), where 6 is the initial angle 
of magnetization whose distribution is given by Eq. (1). The corresponding distribution of 
switching time is calculated from the definition: D(t)cIt = —P{9) smOdd. By defining the 
switching probability density D{r)dr = — P(n m , 0) sin8d8 and utilizing the above relation 
between r and 6, the probability of not being switched is 

Jo 4 V. r 

where a long time limit is taken. The relaxation rate is identical with Eq. (17). This suggests 
the initial-condition randomization is the leading cause for switching time distribution in 
the limit of large £ and I/I c ^> 1. Deviation occurs when £ is small, such as shown in Fig. 1, 
or when the current is near I c , as in Fig. 2. 
This work is supported in part by IBM. 
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Figure Caption 

FIG. 1 (Color online) Numerical results of the smallest nonzero eigenvalues of Eq. (16) for 
increasing £, in cases of different I /I c . Tq 1 = a^(H + H k ). 

FIG.2(Color online) Numerical results of the smallest nonzero eigenvalues of Eq. (16) for 
increasing I / I c — 1, where £ = 80. The red line is the fitting line: Air = 2(7 / I c — 1). 
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